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Abstract. Millimeter interferometry provides evidence for the presence of mm to cm size 'pebbles' in the outer 
parts of disks around pre-main-sequence stars. The observations suggest that large grains are produced relatively 
early in disk evolution (< 1 Myr) and remain at large radii for longer periods of time (5 to 10 Myr). Simple 
theoretical estimates of the radial drift time of solid particles, however, imply that they would drift inward over 
a time scale of less than 0.1 Myr. In this paper, we address this conflict between theory and observation, using 
more detailed theoretical models, including the effects of sedimentation, collective drag forces and turbulent 
viscosity. We find that, although these effects slow down the radial drift of the dust particles, this reduction is not 
sufficient to explain the observationally determined long survival time of mm/cm-sized grains in protoplanetary 
disks. However, if for some reason the gas to dust ratio in the disk is reduced by at least a factor of 20 from the 
canonical value of 100 (for instance through photoevaporation of the gas), then the radial drift time scales become 
sufficiently large to be in agreement with observations. 
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1. Introduction 

Millimeter and sub-millimeter observations have shown 
the presence of large amounts of millimeter to centimeter- 
sized grains in the outer regions (~ 100 AU) o f disks 



around Herbig Ae and T Tauri star s iTestiet al. . 2005 : 



Rodmann et all . I200d IWilner et all . l2005l : iNatta et al 
2007^ . The presence of these grains, which are much 
larger than the grains typically found in the interstel- 
lar medium, is often regarded as evidence that the first 
steps of planet formation are taking place in these disks. 
The presence of such large grains, however, also poses 
a major problem. According to simple theoretical con- 
sider ations, grains of this size undergo a rapid radial 
drift ( Whippl j . 1972 : Weidenschilline . 1977 ). which causes 
them to disappear from the outer disk in a very short 
time (iTakeuchi fc Linl. l2005l: iKlahr fc Bodenheimeii l2006t 



Alexander fc Armitagei r 2007 ). However the typical age 
of protoplanetary disks that are observed at millimeter 
wavelengths is a few Million years, which is much longer 
than this radial drift time scale. Takeuchi & Lin pro- 
pose that either the entire grain growth process is slow 
or that the grains are the collision products of a popula- 
tion of even larger bodies (> 10 m). The second expla- 
nation requires that in addition to the grain population 



that is observed at mm wavelengths, there is a population 
of smaller/larger bodies which act as a reservoir of solid 
material from which mm/cm-sized grains are continuously 
produced. The problem is that if the drift time scale is, 
for example, 20 times shorter than the disk life time, this 
reservoir of larger bodies must contain at least 20 times 
more mass than the observed dust mass. If it is assumed 
that the particle size distribution follows a powerlaw then 
the total mass of the disk and the minimum upper par- 
ticle size of this distribution can directly be calculated 
from the slope of the mm flux of the protoplanetary disk. 
This analysis shows that the amount of dust responsible 
for the millimeter fluxes of these disks is in many cases 
already very high^ of the order of 10~^ Mp, or even higher 
( Testi et all 120 03: Nat ta et al.l. 120041: IWilner et all . 12005: 
Rodmann et al.. . ,2006 : iRodmanni " 20061 ). A 20 times more 
massive reservoir of larger (non-observable) bodies is then 
clearly unrealistic. These arguments suggest that perhaps 
the standard theoretical estimate of the radial drift may 
be not applicable. 

Interestingly, a related drift problem is present for the 
theory of planetesimal formation. The drift problem for 
mm/cm-sized particles at 100 AU, that will be investi- 
gated in this paper, similarly shows up for meter-sized 
particles at 1 AU. The radial drift of these large bodies 
in the inner parts of the disk is so high that they should 
drift into the evaporation zone over time scales of ^ 10^ 
yrs. This is one of the fundamental and unresolved prob- 
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lems of planet formation ( Dominik et al. . 20071 ). In that 
sense, the cm problem at 100 AU is a proxy for the meter 
problem at 1 AU, and figuring out a solution at 100 AU 
may give important clues to what happens at 1 AU. 

The goal of this paper is to study the radial drift of 
mm/cm-sized particles in more detail. We will investi- 
gate the magnitude of the drift problem, and which effects 
might keep the grains for a few Myrs in the outer parts 
of the disk. To address this issue we will proceed in steps. 
In the first st ep we will review the radial drift of in divid- 
ual particles JWhipp ld. 119721 : IWeidenschilhngl . Il977l ) . This 
section will show that the drift time scale of such particles 
is orders of magnitude smaller than the age of the disks 
observed (5 to 10 Myrs). In a second step, we explore 
the possibility that collective effects of the dust might 
slow down the drift. Collective effects tak e place when the 
dust settles into a thin m i dplane layer dOubrulle et al 



1995t iGaraud et al] . 12004 ISchrapler fc Hennind . 12004 ) 
This process increases locally the dust-to-gas ratio, and 
the dynamics of the dust starts to affect the gas motion 
( Nakagawa et al.l . ll986l : IJohansen et"al] . l2006al ). This may, 
in turn, reduce the relative velocities between the dust and 
the gas, and hence reduce the head wind that causes radial 
drift. We will investigate the magnitude of the reduction 
and if thin midplane layers yield a possibility to increase 
the drift time scales to some Myrs. In the third step we 
improve on these calculations by including vertical angu- 
lar momentum exchanges in the disk through turbulent 
viscosity. Finally, we speculate on other potential ways in 
which mm/cm-sized grains could be prevented from drift- 
ing inward on a time scale shorter than the life time of the 
disk: particle trappi ng in yortices and gas pressur e max 



ima 



(iBarge 



19951: iKlahr fc Hennind . Il997 



Froma,ng fc Nelsonl. l2005t iJohansen et al.l . l2006bl ). spiral 



arms ( Rice et al. . 2004 ) and photoevapor ation of the gas 



leaving the dust behind (jAlexander et al K2OO6) 

Over the last decades various papers determined the 
radial drift of dust particles and the structure and dy- 
namics of thin midplane dust layers. In particular the lat- 
ter problem has attracted much attention, but for an en- 
tirely different reason than ours: Gravitational instabili- 
ties in thin midplane dust laye rs are thought to b e a pos - 
sible origin of planetesimals (jGoldreich fc Wardl . Il973f ). 
A lively debate has since appeared about the viability 
of this concept, spurring various papers including mod 
els of midplane dust lay e rs (IWeidenschil ling. 1980; S ekiva 



1998; 
2006; 



Weidenschilline 



Youdin fc Sh 



■20061 : ICuzzi fc Weidenschilhne . 
2OO2I : lYoudin fc Chiand . 12004 : 



Johansen et al.l . l2006al) . The richness of this literature 
gives an indication of the complexity of the problem. 
Hence, due to this complexity only a sub-set of the possible 
physical effects are considered in these dust layer models. 
In particular the collective effects of the dust and the ef- 
fects of vertical and radial viscosity have not been studied 
yet in combination. Therefore, another goal of this paper 
is to present a model of dense dust midplane layers that 
include a multitude of physical effects, albeit still in the 
form of a 1-D vertical model. 



In this paper we do not consider dust particle co- 
agulation. Larger particles in protostellar disks beyond 
sizes that can be found in the interstellar medium form 
by collisional sticking due to relative motion s of the dust 
(|Beckwith eralll2000l : lBlum fc Wur^ .l200rf). However, at 
100 AU in the disk the dust particles rapidly drift inwards 
before they can even grow to sizes that are discussed this 
paper (Brauer et al. in prep.). In the present paper we 
will ignore the issue of the formation of these mm/cm 
size grains and assmne that all the dust is in the form of 
grains of a given size (a 'monodisperse' size distribution). 
We will then investigate whether these grains can remain 
in the outer regions of the disk for a sufficiently long time. 

The structure of this paper is as follows. In section[2]we 
discuss the radial drift of the dust and its radial drift time 
scales. Section [3] includes other possibilities to increase the 
drift time scales. A detailed discussion of the theoretical 
background is given in the Appendix A. 

2. Radial drift of the dust 

We consider a T Tauri star with a stellar mass of 1 Mq, 
a stellar radius of 2.5 Rq and a surface temperature of 
4000 K. We assume that the disk around this pre-main 
sequence star h as an inner and an ou ter radius of 0.03 
AU and 150 AU (|Rodmann et al.Ll2006h . respectively. The 
mass of the disk Mdisk is a free parameter of our model. 
We assume that the disk is passive and irradiated by the 
central star under an angle of 0.05 rad. At a disk radius of 
1 AU these values imply a midplane temperature of 200 
K assuming that the disk is isothermal in the vertical di- 
rection. Moreover, we will consider a turbulent disk. The 
amount of turbulenc e in the disk is descr i bed b y the tur- 
bulent a parameter ( Shakura fc Sunvaev . 19731 ). which is 
a free parameter of our model. 

Now, we will calculate the radial drift of the dust and 
the radial drift time scales. In order to demonstrate the 
physics behind the calculations and its implications on the 
radial velocities we will proceed in certain steps. In every 
step more effects are included to demonstrate the influence 
on the drift velocities (cf. Fig.[T]). 

1. The first subsection addresses the drift of individual 
particles in a gaseous disk. These results are valid when 
the dust and the gas are well mixed. 

2. However, under certain conditions (low-turbulence 
disks, large dust particles), the dust sediments so close 
to the midplane that the density of the dust exceeds 
the density of the gas. In that case collective effects 
of the dust come into play, which is the topic of the 
second subsection. 

3. Finally, in the last subsection, we also include the effect 
of vertical angular momentum exchange between the 
dust midplane layer and the gaseous layers above the 
midplane. 

We are primarily interested in drift time scales needed to 
keep mm to cm size particles in the outer parts of the 
disk, hence, we will calculate characteristic drift times at 
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Fig. 1. This figure shows the three scenarios we will inves- 
tigate in this paper. In step 1 we consider the drift of in- 
dividual particles. In step 2 we investigate the influence of 
collective effects on the particle drift. These effects become 
of importance if the particles settle into a thin midplane 
layer. Finally, in step 3 we will investigate what happens 
if we include turbulent viscosity in our simulations. In this 
case, angular momentum is exchanged between the dusty 
midplane layer and the gaseous layers above the midplane. 
The horizontal arrows indicate the radial velocity of gas 
and dust. The curved arrows indicate the azimuthal ve- 
locities. 



the very end of every of these three subsections. Wc will 
see how these time scales are affected by including more 
and more effects. 



2.1. Step 1 - Radial drift of individual particles 
2.1.1. Equations 

The fundamental cause for inward drift of the dust is the 
difference in velocity between gas and dust. While the dust 
moves with Keplerian velocity the gas moves slightly sub- 
Keplerian. This is due to the fact that the gas is not only 
affected by the gravitational and the centrifugal force but 
additionally feels a radial pressure force that does not act 
on dust particles. This extra force is caused by the de- 
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Fig. 2. The Stokes number St as a function of particle 
radius a at different radii in the disk. In this calculation 
the solid particle density of the dust is 1.6 g/av? and the 
stellar mass is 1 Mq. The innner and the outer radius of 
the disk are 0.03 AU and 150 AU, respectively. The disk 
mass is lQ~'^Mp,. 



crease of gas pressure in the radial direction. Since this 
force, which exclusively acts on the gas, partly compen- 
sates gravitation, the gas moves slower than Kepler speed 
and therefore slower than any dust particle in the disk. 
Hence, the dust particle feels a continuous headwind from 
the gas. This headwind causes the dust particle to lose its 
angular momentum and spiral inward. 

Whipple! (|l972h formulated the first equations for 



the radial drift of ver y small and very large particles. 
Weidenschillind (1977) later derived a set of equations 
with a general drag force to calculate the radial drift of 
solid particles of any size. We will formulate all equa- 
tions in the dimensionless Stokes number formulation. The 
Stokes number can be regarded as a measure of grain size 
(see Fig. [2]) . The full definition is given in Appendix IA.4I 
In terms of this dimensionless formulation, these equations 
aquire the forrr0 



St 

,,2 



T 



wz. 



(1) 



The variables Wr and w^p denote the radial and azimuthal 
velocity of the dust, respectively. The form of the drag law 
is implicitly included in the Stokes number. The quantity 
wn is the velocity by which the gas moves azimuthally 
slower than Keplerian velocity Vk, i.e. Wgas = Vk — f n- The 
velocity vtsi will also turn out to be the maximum radial 
drift velocity of the dust. We will take wn from here on as 
our " standard velocity" scale (cf. Eq. IA.22P apart from the 
Keplerian velocity V\^. For our disk model the quantity v-^ 



^ The physical meaning of all variables is listed at the very 
end of the paper. 
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Fig. 3. The radial drift velocity of individual particles Wr 
in units of vn as a function of Stokes number. 



is given by ( Weidenschilling . 1977 : Nakagawa et al. . 1986h 



= 1.28 



(2) 



In this expression Cg denotes the isothermal sound speed. 
The number 1.28 is due to the specifications of our disk 
model (cf. App. IA.6|) . Note that the velocity scale wn does 
not depend on the disk mass. In our disk model this quan- 
tity is also not dependent on the location r in the disk. 

The set of equations ^ is generally difficult to solve 
and only numerical methods provide information about 
the drift velocity. However, in some cases the situation 
simplifies. If the Stokes number is not dependent on the 
particle velocity, the equations ([T]) can be solved analyti- 
cally. Assuming this independency a straightforward cal- 
culation yields 

9 

(3) 



V-N ■ 



This equation directly shows that the drift velocity has a 
maximum when the Stokes number is unity and the max- 
imal drift velocity is wn- The drift velocity as a function 
of Stokes number is shown in Fig. 

2.1.2. Drift time scales for individual particles 

We will now investigate the radial drift times of individual 
particles in the outer parts of the disk. More specifically 
we are interested in the conditions on particle radius and 
particle porosity that provide time scales larger than a few 
Myrs. 

The drift timescale Tdrift equals 



Tdrift = , 

Wr 



(4) 



which should correspond to the age rage of the disks ob- 
served, thus a few Myrs. Throughout the paper we focus 
on the radius of r = 100 AU. To find the critical Stokes 
numbers for which the radial drift time scale equals the 
age of the disk, we replace Tdrift in Eq- Q by Tago . We then 



insert Eq. ^ into Eq. and solve for St. This yields 
two critical Stokes numbers: 



Std 



cm 



± 



(5) 



The interpretation of these two numbers is the following: 
If the Stokes number of the dust particle falls into the in- 
terval [St_,St+], then the drift timescale is shorter than 
Tage . If it falls outside of this interval, then the drift time 
scale is long enough that these particles can be observed in 
the protoplanetary disk of age Tago . Since the Stokes num- 
ber interval is a rather abstract depiction we reformulated 
it into a particle radius interval with a similar meaning. 

The region of too short time scales at 100 AU is shown 
in Fig. m as a function of disk mass Mdisk and surface den- 
sity S. In this diagram we applied a dust material densitjH 
Ps — 1.6 g/cm-^, a Kepler frequency fJk = lO^^^/s and 
Cs = 2.6 X 10^ cm/s (corresponding to a temperature of 
20 K) . We take as the age of the disk Tagc = 2 Myrs. The 
maximum radial drift velocity at this location is v-^^ — 60 
m/s. We will make use of these values at all times in this 
paper unless otherwise noted. The two Stokes numbers 
St± that are imphed by these values are St_ — 0.002 and 
St+ — 474, corresponding to the lower and upper edge of 
the grey zone in Fig. [J] respectively. 

The figure shows that the particle radius interval in 
which the time scale of individual particles is shorter than 
2 Myrs ranges over more than 5 orders of magnitude in 
radius. Particles ranging from mm to cm in size are com- 
pletely included in this region independent of disk mass. 
The drift time scale of sub- millimeter particles may exceed 
Tagc when the disk mass is higher than 0.2 Af^. 

The Stokes number as the crucial value for radial drift 
is not only affected by particle radius but a lso by particle 



prope rties like porosity or fractal growth (jKempf et al 
20001 ). This effect of noncompact growth may be consid- 



ered by introducing the filling factor of the particle / de- 
fined by Trip = VpPsf, where mp and Vp are the mass and 
the volume of the particle, respectively. In Fig. (U we also 
calculated the critical particle radius interval for a filling 
factor of f ~ 10^^ (dotted lines). 

The lower filling factor shifts the critical particle radius 
interval towards higher particle radii. The drift time scale 
of mm size particles exceeds 2 Myrs when disk masses 
higher then 0.2 M^, are considered. For cm size particles 
the time scale never exceeds 2 Myrs. For filling factors 
lower than 10~^ the drift time scales of mm and cm size 
particles exceed 2 Myrs for any disk mass higher than 
10^'^ Mi,. However, particles of mass 1 g and filling factors 
of 10^"^ would imply a particle diameter of 5 cm. Since 
this particle size f alls into the regime wh ere compaction is 
thought to occur (jBlum fc Wurm . I2OOOII this filling factor 
represents an unlikely case. 



^ 10% silicate, 30% carbonaceous material and 60% ice 
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following set of equations ( Nakagawa et al. . 1986f ) 
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Fig. 4. The particle radius interval in which the individual 
drift time scale at r —100 AU becomes shorter than 2 
Myrs as a function of surface density or equivalently disk 
mass (shaded region) . The disk mass is computed from the 
surface density using Eq. (|A.21|) . To illustrate the effect 
of noncompact growth we also calculated the interval for 
a fiUing factor of 10^^ (dotted lines). 



2.2. Step 2 - Collective effects 

The scope of the previous subsection can be expanded by 
including the back-reactions from the dust to the gas. We 
do no longer consider a single particle, but include how the 
entire swarm of dust particles can affect the gas motion. 
The modified gas motion has the effect of reducing the rate 
by which the gas extracts angular momentum from the 
dust, and thereby reduces the radial drift of the dust. Such 
collective effects play the strongest role for low a-values 
so that a thin midplane dust layer can form in which the 
dust density is high. This scenario of reducti on of radial 
drift was described bv iNakagawa et aL ( 19861 ). 



The necessity to take this additional effect into ac- 
count may be illustrated by regarding the following ex- 
treme scenario. We consider a hypothetical disk in which 
the dust density is much higher than the gas density. In 
such a dust-dominated disk the dust is hardly influenced 
by the gas. The gas, which tends to move sub-Keplerian, 
is dragged along with the dust since it feels a continuous 
tailwind of dust particles. Therefore, the gas perpetually 
gains angular momentum from the dust and spirals out- 
ward. The radial drift of the dust is negligible because the 
dust-to-gas ratio is much higher than 1 . This is the reverse 
situation of the case described in Subsection 12.11 In gen- 
eral, though, we neither have a perfectly dust-dominated 
nor gas-dominated situation. We then have to solve the 





— Apd[Ur - 


^ 1 .1 

- Wr ) OrP 

Pg 




- Apd{u^ 


- W'v) 




— Apg{Wr 


- Ur) 




- Apg{w^ 





(6) 



The quantities u and w denote the velocity of the gas and 
the dust in a Keplerian comoving frame, respectively. The 
subscripts r and (p indicate the radial and the azimuthal 
components of the velocities. The variables pg^d denote 
the mass densities of gas and dust and the quantity A is 
defined as A = Jlk/PgSt. We are primarily interested in 
the radial dust velocity of the Nakagawa-Sekiya-Hayashi 
solution (NSHs) which can be expressed as 



,NSH 



and 



1 



1 



(7) 



The quantity e = pd/pg denotes the local dust-to-gas ratio. 
When ip ^ 1, i.e., when the dust-to-gas ratio is approach- 
ing zero, Eq. ([7]) reduces to the corresponding equation 
for single-particle drift, Eq. ([3]). The difference between 
these two equations is the additional factor ^ that modi- 
fies the Stokes number St and the maximum drift velocity 
wn. Since ijj is always smaller than one, the collective radial 
drift of the dust will always be smaller than the individual 
particle drift. 

Taking collective effects into account requires knowl- 
edge about the dust density. Therefore, certain disk pa- 
rameters, i.e. the turbulence parameter a as well as the 
initial dust-to-gas ratio eo before sedimentation, become 
of importance. Another dimensionless number, which now 
comes into play, is the turbulence parameter q. This num- 
ber determines whether turbulent diffusion is realised 
by small turbulent eddies moving fast or by big eddies 
moving slow (cf. Appendix IA.1|) . These quantities deter- 
mine the thickness of the dust layer h and, hence, the 
dust density (cf. Eqs. IA.16I and IA.18|1 . In these equa- 
tions we assume that the dust density in the verti- 
cal direction has a gaussian shape. This ansatz might 
be put into question if the tur bulence is self- i nduced 
(| Weidenschillingl . 1 1979h . Although ljohansen et all (|2006ah 
showed that the vertical dust density in self-induced tur- 
bulence has a gaussian shape for canonical dust-to-gas 
ratios, the vertical structure can show a different shape 
especially w hen e is increased for in stance through photo- 
evaporation ( Weidenschillingl . 20061 ). 

Since the dust density is a function of height above the 
midplane z, the radial drift velocities are dependent on z 
as well. This vertical dependency is shown in Fig. ([5]) for 
an exemplary NSH solution. In this calculation we applied 
the values St = 1, a = 10-^ q = 1/2 and eo = 10"^ 

The plot shows that the dust moves inwards while the 
gas moves outwards which is generally the case in the NSH 
solution. In the higher regions of the disk {\z\ > 3 h) the 
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Fig. 5. The upper figure shows the collective radial veloc- 
ities of gas and dust of the laminar NSH solution in terms 
of wn as a function of height above the midplane. The 
lower figure shows the radial mass flux of gas and dust in 
arbitrary units. The values applied in this calculation are 
St 1, a = IQ-^, q = 1/2 and eq = 10"^ 



dust-to-gas ratio is much smaller than unity causing the 
collective drift behaviour to match the individual particle 
drift. However, closer to the midplane of the disk collective 
effects become important. With increasing dust-to-gas ra- 
tio towards the midplane, the radial inward drift of the 
dust decreases while the gas starts to move outwards. The 
clear difference in velocities between the collective drift 
and the individual drift around the midplane along with 
the fact that most of the dust is located in this region 
demonstrates the importance of collective effects for disks 
with low turbulence. 

The radial velocities as a function of height above the 
midplane do not directly tell something about the entire 
radial flow of the dust since the dust itself is vertically dis- 
tributed in a certain way. For this reason we will now cal- 
culate the vertically averaged radial velocity of the dust. 
This integrated velocity is given by 



-NSH 



(8) 
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A contour plot of this quantity as a function of the tur- 
bulence parameter a and the Stokes number St is shown 



10"= 10"' lO"-" 10"' 

Turbulence parameter a 

Fig. 6. A contour plot of the vertically averaged collective 
radial drift velocities of the dust in terms of the individual 
particle drift velocity as a function of turbulence parame- 
ter a and Stokes number St. The numbers in the diagram 
indicate the contour level and are related to the line on 
the left, respectively. The parameters for this calculation 
are q = 1/2 and eo = 10"^. 



in Fig. The drift velocities in this diagram were ex- 
pressed in terms of the corresponding individual particle 
drift velocity in order to explicitly point out the differences 
between these two models. 

The figure shows that for fixed Stokes numbers the de- 
viation increases continuously with lower turbulence in the 
disk. Lower a values imply thinner dust layers and, there- 
fore, higher dust-to-gas ratios. With higher dust-to-gas 
ratios the back reaction of the dust to the gas increases, 
and hence the deviation between individual and collective 
drift velocities. One obvious solution to the whole radial 
drift problem of grains in the outer parts of the disk would 
be to continuously decrease the amount of turbulence in 
the di sk or even to set a to zero. However. I Weidenschillinel 
' las shown that a shear- instability between the dust 
layer and the gas induces a weak, but non-negligible level 
of turbulence. This turbulence is called 'self-induced tur- 
bulence' which constrains the a- value to be at least of the 
order of ~ 10~^ (cf. Appendix [B|) . 

For fixed a Fig. [S] shows that for low Stokes numbers 
(small grains) the drift behavior approaches the individual 
particle drift. Low Stokes numbers imply thick dust lay- 
ers, causing low dust-to-gas ratios. For high Stokes num- 
bers (large grains), very thin dust layers are obtained. One 
would intuitively think that this maximizes the collective 
effects. However, as can be seen from Eq. ((T)), in the limit 
of St oo one gets u^^^ 2v-f^/St which is equal to 
the individual particle drift of Eq. ^ . So for large St the 
radial drift indeed drops, but not due to collective effects. 

For a Stokes number of unity and a turbulence param- 
eter of 10^® the dust-to-gas ratio in the midplane is given 

by 

Cmid = £o-^ = eo V — = 10. (9) 
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Fig. 7. As Fig. 21 but now with collective effects of the 
dust and the gas included, at r = 100 AU. Different grey 
scales are the results for different levels of turbulence, 
hence different thicknesses of the dust midplane layer. 
Note that for a = 10^^ the results are virtually identi- 
cal to the single-particle case (Fig. S]), since in this case 
the dust layer is so thick that the dust-to-gas ratio is much 
less than unity. 



Therefore, the collective radial drift in the midplane in 
terms of the individual particle drift according to Eq. ^ 
is wJ^fj^/uN = (1 + e)"^ = 0.008. However, the vertically 
averaged drift velocity of the dust in terms of the individ- 
ual particle drift in Fig. ^ is only 0.07, which is almost 
one order of magnitude higher. The reason for this is that 
the largest radial dust mass flux is not in the midplane, 
but slightly above the midplane (see Fig. (5]) . The mass flux 
is the product of dust density pd and dust radial velocity 
yNSH^ Although Pd drops strongly slightly above the mid- 
plane, the radial velocity u^^^ increases even faster, so 
that the product pd u^^^ has a maximum slightly above 
the midplane. While the formation of a dense dust layer 
can reduce the radial drift velocity in the midplane by a 
factor of 100, the vertically averaged radial drift velocity 
can be only reduced by a factor of at most 10. 



2.2.1. Radial drift times including collective effects 

Armed with the above drift velocity expressions we now 
calculate the conditions on particle radius and particle 
porosity that provide time scales larger than 2 Myrs taking 
into account collective effects. 

At first we focus on the conditions on the particle 
radius. The interval of this particle property that corre- 
sponds to time scales shorter than 2 Myrs is shown in 
Fig. (O for different a-parameters. The second turbulence 
parameter q is fixed at 1/2 at all times, the initial dust- 
to-gas ratio is 10~^ and the filling factor / is unity. All 
other parameters were already mentioned in the last sec- 
tion and are not changed throughout the paper unless di- 
rectly stated. 




Mdisk [ M. ] 
0.010 0.100 



0.100 



0.010 



0.001 r 



t < 2 Myrs 




t > 2 Myrs 



Z [ g/cm^ ] 

Fig. 8. The dust particle filling factor that provides time 
scales larger than 2 Myrs as a function of disk mass for 
different turbulence parameters at r = 100 AU in the disk. 
In this calculation collective effects of dust and gas are 
taken into account. 



According to this plot, the critical particle radii that 
provide the requested time scales hardly differ from the 
critical particle radii of the individual particle drift calcu- 
lated in the last section. Even for small turbulence param- 
eters which favour collective effects the time scales for mm 
to cm size particles are shorter than 2 Myrs for any disk 
mass considered. The reason for this is that for very high 
and very low Stokes numbers, like the two critical numbers 
St_ and St+ representing the boundaries of the grey areas 
in Fig[71 collective effects play a minor role (see Fig.El and 
discussion in the last subsection). The Stokes numbers for 
which the collective effects play the strongest role lie in 
the middle of these grey areas, i.e. where the drift time 
scales are anyway much too short to be compatible with 
the observations of mm-sized particles in protoplanetary 
disks. 

So what about the effects of fractal or porous growth? 
For simplicity we set the mass of the dust particle to be 
1 g and then calculate the particle filling factor that pro- 
vides time scales larger than 2 Myrs. For a filling factor 
of unity a dust particle of 1 g corresponds to a particle 
radius of 1/2 cm. For / = 10^"* the particle radius can 
be calculated to be 11 cm. Like in the last paragraph we 
will perform the time scale calculation in dependency on 
the disk mass. The results of this calculation is shown in 
Fig. ([5]) for different values of the turbulence parameter 
a. This diagram shows that for filling factors lower than 
10~^ the time scale exceeds 2 Myrs subject to the con- 
dition that the disk mass is higher than ~ 0.2 M^,. For 
even higher disk masses the filling factor may exceed 0.1 
for certain turbulent a parameters. T his filling factor cor - 
responds to a particle radius of 1 cm. lOrmel et al.l (j2007l ) 
showed that particle growth in protostellar disks can be as- 
sociated with filling factors of less than ^ 10^^. Therefore 
fractal growth seems to be an actual possibility to consid- 
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erably increase the radial drift time scales of the dust. We 
will come back to that point in the discussion. 

2.3. Step 3 - Effect of turbulent viscosity 

We will now investigate the role of turbulent viscosity in 
addition to the effects studied so far. Including viscosity 
terms will have the opposite effect on the drift velocities 
than the collective effects. It will increase the radial drift 
of the dust and shorten the drift time scales. We will give 
the Navier-Stokes equation (NSe) including collective ef- 
fects and viscosity terms and solve these equations numeri- 
cally. However, we would like to discuss first why turbulent 
vicosity increases the radial drift of the dust. Under cer- 
tain conditions, i.e. small turbulence parameters a or high 
Stokes numbers St, the previous sections have shown that 
the dust may settle into a thin midplane layer. When the 
dust-to-gas ratio inside this layer exceeds unity the gas is 
dragged along with the dust. Both components, dust and 
gas, tend to move with Keplerian velocity. 

On the other hand, above the dust layer the dust-to- 
gas ratio is much smaller than unity. In this region the dust 
particles feel a continuous head wind which forces them 
to move with slightly sub-Keplerian velocity. This vertical 
decrease in azimuthal velocity from Keplerian velocity in 
the midplane to sub-Keplerian velocity in the higher re- 
gions of the disk produces a nonlinear velocity gradient in 
both the gas and the dust. 

If we include viscosity in our considerations it attempts 
to damp nonlinear spatial velocity differences. The vertical 
velocity gradient described above represents such a differ- 
ence. Turbulent viscosity now acts in such a manner that 
it transports angular momentum from the midplane to the 
higher regions of the disk. While the midplane, the region 
where most of the dust is located, loses angular momen- 
tum and falls inward, the regions above the midplane gain 
angular momentum and move outward (cf. Fig. [T]). This 
mechanism of verti cal angular momentum e xchange was 
first investigated by lYoudin &: Chiang! (|2004l ). 



2.3.1. Navier-Stokes equations 

The Navier-Stokes equations for this problem are basically 
the set of equations ^ plus some second order derivative 
terms due to the inclusion of viscosity 



= 




— Apd{Ur — Wr) 


= 




~ Apd{u^ ~ w^) 


= 




— Apg{Wr — Ur) 


= 


-ifikWr 


- Apg{Wip - U^) 



1 



drPg 



— dz (pdOzWr) 
Pd 

—dz (pddzW^) . (10) 

Pd 

Hence, the algebraic Eqs. ^ turn into four coupled, dif- 
ferential equations of second order. The left hand side of 
the Navier-Stokes equations representing the time depen- 
dencies are set to be zero since we are interested in steady 



state solutions. The vectors u and w denote the velocities 
of the gas and the dust, respectively. The first terms on 
the right side correspond to the Coriolis force. These terms 
arise from the fact that the equations are formulated in 
a comoving frame. The second term represents the drag 
force coupling between the gas and the dust. The effects 
of viscosity show up in the third terms. The expressions 
for the viscosity of the gas and the dust can be found 
in Appendix [X] In the following we will denote viscos- 
ity terms which involve derivatives of radial (azimuthal) 
velocities as 'radial (azimuthal) viscosity terms'. The very 
last term in the first line is an extra force acting on the gas 
which is caused by a radial pressure gradient. This term 
is responsible for the gas moving slower than the dust and 
causes the radial drift. The densities of gas and dust serve 
as input for the Navier-Stokes equations. 



Takeuchi &: LinI (j2002i ) also investigated the effect of 
gas viscosity on the drift of dust particles, but they ne- 
glected collective effects. This allowed to solve the equa- 
tions analytically. The drift of the dust particles in their 
calculations was a superposition of two different effects: 
The individual dust particle velocity with respect to the 
gas and the velocity of the gas itself. 

The former part of the dust particle drift was discussed 
in detail in Section ITTl The second part of the dust parti- 
cle drift investigated by Takeuchi and Lin was due to the 
gas accretion process. This process of the gas is associated 
with a certain radial accretion velocity. Since the dust is 
to some extent coupled to the motions of the gas the dust 
is carried along with the accreting gas which leads to an 
extra source of radial particle drift. 

In the beginning of this section we described that gas 
viscosity also increases the radial drift of the dust when 
collective effects come into play, i.e. when the dust settles 
into a thin midplane layer and starts to affect the motions 
of the gas. This process is different from the single particle 
considerations discussed by Takeuchi and Lin since it is 
caused by collective effects and not by gas accretion. In 
the following, we will estimate the ratio of these two radial 
drift velocities. 



The additional drift due to the accretion process may 
be estimated by a characteristic acc retion velocity of the 
gas w hich is given by Uacc oc ac^/Vi^ (jShakura fc Sunvaev . 
19731 ). Viscous collective effects imply dr i ft vel ocities of 
order Wcoii = c^/V^eoRe'' (| Weidenschillind . l2003l ) . The ra- 
tio 



(11) 



has values of at most « 10 ^ foJl a = 10 ^, eq = 10 ^, 
Re* = 10^, Vk = 3 X 10^ cm/s and Cs = 3 x 10^ cm/s. 
For smaller turbulence parameters which we will consider 
in this paper the influence of gas accretion will be even 
lower. 
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Fig. 9. The azimuthal and radial gas and dust velocities 
as a function of height above the midplane. The dotted 
lines denote the analytical normalised NSH solution with- 
out viscosity {ly = 0) as discussed in Section The solid 
lines indicate the numerical solution of the Navier-Stokes 
equations including azimuthal and radial viscosity terms. 
The dashed line in the top diagram shows the radial veloc- 
ity of the gas if radial viscosity terms are neglected. The 
values for this simulation are St = 1.2, a = 10~^, q = 0.5 
and eo = 0.01. Note that more than 80% of the dust in 
within the z— interval [~2h,2h]. 



2.3.2. Numerical results 

The parameters for the simulation are St = 1.2, a = 10^^, 
g = 1/2, eo = 10"^. The results are shown in Fig. O The 
two parameters a and St are chosen in a way that effects 
of viscosity become visible. These values imply a dust- 
to-gas ratio of 5 in the midplane of the disk and a half 
thickness of the dusty midplane layer which is w 0.002 H. 
The turbulent motions of the gas have a speed of 0.005 
UN- The dotted lines in the Fig. [5] indicate the analyti- 
cal solution of the laminar NSH equations {i^ — 0) which 
were discussed in the last section. The solid lines in Fig.IH 
indicate the numerical solution including viscosity which 
differs significantly from the NSH solution. 

Let us focus on the radial dust velocity since we are 
primarily interested in radial drift time scales. The radial 
flow of the dust is significantly affected by turbulent vis- 
cosity if a is smaller than 10^"*. In this regime the effect 
is largest for Stokes numbers w 5. For a > 10~* the radial 
flow is approximately the flow predicted by Nakagawa et 
al. and viscosity seems to play a minor role. The radial 
dust velocities in the midplane with and without viscosity 
terms may differ by a factor of 5 for small a parameters 
and St w 5. 

The azimuthal dust and gas velocities as a function 
of height above the midplane vary in a complex manner. 
However, for Stokes numbers smaller than unity the sit- 
uation with regard to the azimuthal velocities simplifies. 
In this regime these velocities do hardly differ from the 
expression given by Nakagawa et al. and viscosity seems 
to be negligible. 

The radial outflow of the gas, which is shown in 
Fig. is reduced if turbulent viscosity is included. This 
decrease may be up to a factor of 30 for small St and a 
parameters. For turbulence parameters higher than 10^^ 
the radial net outflow of the gas differs less than 10% from 
the outflow predicted by the NSH equations. 

2.3.3. Width of the azimuthal gas velocity layer 

The calculation of radial drift velocities in protostellar 
disks including collective effects and effects of viscosity 
are a challenging topic. Most equations can not be solved 
analytically and only numerical solutions provide informa- 
tion on the evolution of these disks. Therefore, disk model 
simplifications often come into play. 

One simplification is that the dust sub-disk is assumed 
to be extremely thin and thought to behave, to some ex- 
tent, like a solid disk. This approximation i s calle d "plate 
drag approximation" (jGoldreich &: Wardl . Il973[ ). Under 
this condition, the gas layer above the dust layer can 
be described by an Ekman layer: The gas in the mid- 
plane is forced to move along with the Keplerian rotat- 
ing solid equatorial subdisk. High above the midplane the 
gas is in equilibrium with the radial gas pressure gradi- 
ent, yielding a slightly sub-Keplerian rotational velocity. 



Estimated values at 100 AU 
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Fig. 10. The width of the azimuthal gas velocity distri- 
bution in units of the dust scale height /i as a function of 
Stokes number St for 3 different turbulence parameters a. 
The dotted line indicates the ~ St ' dependency of the 
Ekman layer. 



The Ekman layer is the transitional region between these 
two extremes. The thickness of this layer depends on the 
viscosity of the gas. 

In this subsection, we will compare our results with 
the predictions of the simplified model described above. 
We want to know the extent of the region where gas and 
dust affect each other and effects of viscosity become of 
importance. The comparison with regard to the drift ve- 
locities implied by this approximation, however, will be 
discussed in Section [2.3.61 

To quantify the length scale over which viscous collec- 
tive effects play an important role, we define a measure A 

by 

A(a,St) = y g{z)\z\ dz. (12) 

The function g(z) is given by 

g{z) = cn|u^ + vj<i\z'^ . (13) 

The constant cn provides the normalization of g. With this 
distribution function, deviations from the single particle 
solution u^p+v-^ are weighted in a way that differences high 
above the midplane are more important than differences 
close to the midplane. Therefore, A provides informations 
about the width of the vertical azimuthal velocity distri- 
bution of the gas. The dependence of this quantity as a 
function of St is shown in Fig. (fTU|) for 3 different a- values. 

According to this diagram the value of the length scale 
A is a few dust scale heights as long as the Stokes number 
is smaller than unity. For higher St values this quantity in- 
creases exponentially up to more than 40 h for St = 100. 
We also find that A is hardly dependent on the turbu- 
lence parameter a. This behaviour may be understood by 
investigating the length scale of the classical Ekman layer 

^-lxf¥-^T-^^- (14) 




-15 -10 -5 5 10 15 
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Fig. 11. This graph shows dzU^ which indicates the flow 
of angular momentum in the vertical direction. The solid 
and the dotted lines show the flow of momentum with and 
without radial viscosity terms, respectively. The azimuthal 
viscosity terms are included at all times. 



The dotted line in Fig. (fTO|) indicates this dependency 
which shows that the gas layer in fact acts like an Ekman 
layer when the Stokes number exceeds unity. 

2.3.4. Vertical flow of angular momentum 

A remarkable effect of the viscosity is the radial inward 
drift of the gas which is impossible in the laminar NSH 
solution. At certain heights above the midplane the gas 
moves inwards (see Fig. |9]at z = ztSh for example). 

To understand this effect, we provide the basic sce- 
nario. The gas in the midplane of the disk dragged by the 
dust moves azimuthally faster than the gas outside the 
dust layer, which causes a vertical velocity gradient. Since 
viscosity tries to equalize such velocity gradients the gas in 
the higher regions of the disk is accelerated, decelerating 
the gas in the midplane. Therefore, viscosity transports 
angular momentum from the midplane to the higher re- 
gions of the disk. 

To substantiate this effect we calculate the flow of an- 
gular momentum of the gas in the vertical direction. The 
structure of this flow can be analysed by calculating dzU^p 
(see Fig. [TT|) . This calculation was performed with the 
same parameter values as used in the last section. 

The results show that the maximum vertical upward 
flow of angular momentum takes place at ±2/i. It also 
shows that there is a vertical downward flow of angular 
momentum at about ±6/i. This flow is strongly associated 
with the radial inward drift of the gas at certain heights of 
the disk and the inclusion of radial viscosity terms. This 
behaviour can be understood by performing simulations 
without radial viscosity terms: 

If only azimuthal viscosity terms are included the az- 
imuthal gas velocities continuously decrease with increas- 
ing distance from the midplane. This means that angular 
momentum is generally transported in the higher regions 
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Fig. 12. A contour plot of the integrated radial dust ve- 
locity under the influence of viscosity as a function of 
turbulence parameter a and Stokes number St. The num- 
bers related to the contour lines on the left side, respec- 
tively, indicate the net drift velocity in units the individual 
particle drift u. 



Fig. 13. A contour plot of the integrated radial dust ve- 
locity under the influence of viscosity as a function 
of turbulence parameter a and Stokes number St. The 
numbers related to the contour lines on the left side, re- 
spectively, indicate the net drift velocity in units of the 
NSH drift u^^^ without viscosity. 



of the disk and never towards the midplane. This suggests 
that the vertical downward flow of angular momentum is 
an effect caused by radial viscosity. To substantiate this 
assumption Fig. (jlip also shows the vertical flow of angu- 
lar momentum when radial viscosity terms are neglected 
(dotted lines). The inflow vanishes in this case. 

We also calculated the radial velocity of the gas with- 
out radial viscosity terms included. The results of this cal- 
culation indeed demonstrate that the radial inward drift of 
the gas vanishes in this case (see Fig. 1^]). The results also 
show the occurrence of two narrow peaks in the vertical 
velocity distribution of the radial gas velocity without ra- 
dial viscosity terms. These peaks imply high velocity gra- 
dients. Radial viscosity, once included in the simulation, 
reduces these velocity differences by radially accelerating 
the neighbouring regions. This acceleration leads to a de- 
crease in the azimuthal velocities since ^ —u^ due to 
Coriolis forces. This again causes the gas to drift inward. 

Figure [S] also shows that the radial outflow of the gas 
may be faster than if radial viscosity terms are ne- 
glected. The azimuthal velocity differences in gas and dust 
that initially cause any drift behaviour are of the same 
order of magnitude. Therefore, it appears unjustifled to 
neglect radial viscosity terms as often implicitly done by 
using the plate drag approximation for example. 

2.3.5. Integrated radial velocities 

We now calculate the net flow of the dust -S^ according to 
Eq. ([H]). The result is shown in Fig. (fT2|) expressed in terms 
of the individual particle drift velocity. According to these 
results the drift behaviour for high turbulence parameters 
is that of individual particles and neither collective effects 
nor effects of viscosity seem to play a major role in this 
part of the diagram. The net dust velocity has values of 



about —UN for St « 1/2 (cf. Fig. [T^ and decreases with 
lower a values and with growing distance from St « 1/2. 

To demonstrate how viscosity changes the collective 
drift behaviour investigated in step 2, it is suggestive to 
express the viscous collective drift in terms of the NSH 
drift u^^^ . A contour plot of this ratio can be seen in Fig. 
p^ . This plot shows that the radial velocities calculated 
in this section exceed the drift due to collective effects 
by a factor of 2 at most if very low turbulence parame- 
ters are considered. For a parameters higher than 10^* 
viscosity alters the drift scales by a factor of 1.2 in the 
most extreme case. The deviation from individual particle 
velocities due to collective effects were more pronounced 
than those due to viscosity. Therefore, we conclude that 
the drift behaviour is predominantly determined by col- 
lective effects and not by effects of viscosity. 

2.3.6. Plate drag approximation 

Here we will compare our results with previous work. We 
will consider the predictions of the "plate drag approxi- 
mation" . We would like to investigate if these two drift 
models predict the same radial velocities in certain pa- 
rameter regimes. 

In the plate drag approx imation the drift in duced 
by viscosity i s give n by ( Goldreich fc Wardl . 1973 : 
Weidenschillinsl . l2003f) 



CseoR.e*v27r 



(15) 



The derivation of this expression is based on the as- 
sumption that the dust sublayer behaves like a solid disk 
subject to viscous stress on its surface by a turbulent 
boundary layer. This stress extracts angular momentum 
from the dust layer which implies radial drift of the dust 
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Fig. 14. The radial velocity deviation induced by the in- 
clusion of viscosity as a function of Stokes number St for 
different turbulence parameters a. The horizontal lines in- 
dicate the viscous radial drift predicted by the plate drag 
approximation for two different critical Reynolds numbers. 



( WeidenschiUinel . l2006l ). The quantity Re^ denotes the 
critical Reynolds number which indicates the transition 
point between laminar and turbulent flow. The value of 
this number depends on the geometry of the flow and 
is usually determined experimentally. We calculated the 
viscous drift for two different critical Reynolds numbers 
Re* = fOO and Re* = 200. 

We will measure the radial drift due to the inclusion 
of viscosity by the deviation 



-NSH 



(16) 



This difference velocity (5p as well as the predictions of the 
plate drag approximation upd are shown in Fig. (|14p as 
a function of the Stokes number for different turbulence 
parameters a. 

Fig. [m shows that the predictions of these two mod- 
els are roughly of the same order of magnitude if the 
Stokes number is about unity. For Stokes numbers much 
smaller/larger than unity the results of the numerical in- 
tegration of the Navier-Stokes equations deviates from 
the predictions of the plate drag model. Already Youdin 
and Chiang (2004) put the plate drag approximation 
in question. They found that the plate drag overesti- 
mates turbulent stresses and that vertical shear and buoy- 
ancy are important elements missing in this description. 
While the plate drag model involves a radial drift veloc- 
ity which is inversely proportional to the surface density 
of the layer and not e xplicitly dependent on particle size 
Weidenschilih^ (|2006l ) found the very contrary. In his sim- 
ulations, the drift velocity shows no significant variation 
with surface density, but is dependent on particle size 
which clearly speaks against the validity of the plate drag 
model. 



2.3.7. Fitting formula 

A simple fitting formula that reproduces the results might 
be useful for forthcoming purposes for example investiga- 
tions of drift time scales or radial mixing. For this rea- 
son we fitted the vertically averaged radial dust velocities 
given by the numerical solution of Eq. pl?|) . This solution 
includes all effects investigated in this paper, i.e. collec- 
tive effects and effects of turbulent viscosity. This result 
is shown in Fig. (fT5|) . 

For the fit we use an expression of the form 



Ufit 



r(a) 



(17) 



following Eq. ([3]) for individual particle drift. In this ex- 
pression the amplitude F is solely a function of the turbu- 
lence parameter a. The quantity r/2 matches the maxi- 
mum occuring radial dust velocity in units of —I'm when 
a is fixed. The parameter x is given by 



a;(a,St) = lO^^^^Sf^^") . 



(18) 



The two functions ^ and are only dependent on a. The 
parameter ^ determines the location of the maximum of 
the velocity distribution, the parameter p. determines the 
width of the velocity distribution. The fits for the three 
functions F, ^ and /i were performed with polynomials of 
the form 



4 



4 

E 

j=o 



(19) 



in which y is given by y = logj^g Q^- Fo^' the dependency of 
^ on a a second order polynomial turned out to be suffi- 
cient. For the quantities F and /i a fourth order polynomial 
provided a satisfying fit to the simulation results. The co- 
efficients for all these polynomials are listed in table (1). 
The deviation between the fitting function and the simu- 
lation within the parameter intervalls St e [10^^, 10^] and 
a € [10~^, 10~^] is always smaller than 0.01 un- 



i 




c* 







1.89082 


2.06164 X 10-2 


0.57083 


1 


0.14763 


4.69938 X 10"^ 


-0.41644 


2 


0.20912 


-4.05442 X 10-=* 


-0.12910 


3 


8.25120 X 10*^ 




-1.24036 X 10^2 


4 


7.38181 X 10"^ 




-4.09782 X 10"" 



Table 1. Coefficients for the polynomical fit of the simu- 
lated integrated dust velocities. 



2.3.8. Radial drift times including effects of viscosity 

While collective effects reduce the radial drift of the dust 
the additional inclusion of viscosity into the disk model 
again increases it. For this reason, the time scales implied 
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Fig. 15. A contour plot of the integrated radial dust ve- 
locity under the influence of viscosity as a function 
of turbulence parameter a and Stokes number St. The 
numbers related to the contour lines on the left side, re- 
spectively, indicate the net drift velocity in units of — wn- 



by collective effects and effects of viscosity represent an 
intermediate regime between the time scales of individual 
particles and the time scales implied by collective effects. 
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Fig. 16. This figure shows the effect of the dust-to-gas 
ratio on the drift time scales for a turbulence parameter 
of a = 10~^. The figure gives the particle radius interval 
for which the drift time scale is smaller than 2 Myr as 
a function of disk mass and surface density for different 
total vertical dust-to-gas ratios e — Y.d/'Sg- The surface 
density is given at 100 AU. The disk mass and the surface 
density in this figure are the 'original' mass and surface 
density of the disk before the gas depletion that is invoked 
to alter the dust-to-gas ratio. 



3. Other possibilities to increase the drift 
timescale 

We have seen that even with the creation of very dense 
midplane layers for very low a the radial drift is too fast to 
explain the observed millimeter flux of these disks. We will 
now discuss other possible solutions to this problem. We 
will first consider the effect of the dust-to-gas ratio on the 
drift time scales. We will then investigate the importance 
of the turbulence parameter q and consider the possibility 
of non-linear effects which could play an important role. 

3.1. Dust-to-gas ratio 

In this subsection we will investigate the influence of the 
dust-to-gas ratio on the drift time scales. To increase this 
quantity we will remove a certain fraction of the gas from 
the disk. This removal has important implications for the 
drift time scales. When gas is removed from the disk then 
the dust particles are less affected by the motions of the 
gas. This leads to thinner dust layers and hence to higher 
dust-to-gas ratios. For this reason collective effects become 
of importance which reduces the radial drift velocities ac- 
cording to Eq. ([7]). In this paragraph, we will investigate 
how much gas we have to remove from the disk to provide 
time scales larger than 2 Myrs. 

As in the last sections we calculate the particle radius 
interval in which the drift time scale is shorter than 2 
Myr. These calculations were performed for two different 
turbulence parameters a = 10^^ and a = 10~^. The re- 
sults of these simulations are shown in Fig. (fTH|) -(fT7 |l . In 
these figures the disk mass and the surface density are the 



"disk L 1 

0.001 0.010 0.100 




1 10 
Z [ g/cm^ ] 



Fig. 17. This plot is the same as Fig. [HI But here the 
turbulent a parameter is 10~^. 



'original' mass and surface density of the disk before the 
assumed gas depletion that is invoked to alter the dust- 
to-gas ratio. 

The drift behaviour of the dust particles hardly 
changes if only a small fraction of the gas is removed. 
Both figures show that the critical particle radius interval 
is little affected by removing 50% of the gas, cf. Fig. [71 
However, for higher dust-to-gas ratios the critical inter- 
val decreases continuously. Considering the case a = 10^^ 
we find that cm particles are able to stay in the outer 
parts of the disk for 'original' disk masses < 0.05M* and 
> 0.2M^, if only 5% of the gas is left. The critical inter- 
val disappears completely if the total vertical dust-to-gas 
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Fig. 18. This figure shows the effect of the turbulence pa- 
rameter q on the radial drift velocity for two different tur- 
bulence parameters a. The Stokes number was chosen in a 
way that the radial drift velocity for q — 1/2 corresponds 
to a timescale of 1 Myrs with regard to 100 AU assum- 
ing that UN = 60 m/s. For a = 10~^ and a = 10^^ this 
implies St = 225 and St = 235, respectively. 



0.4%. We conclude that q has a very minor effect on the 
radial drift velocities. 

A small h/H ratio leads to a high dust-to-gas ratio in 
the midplane of the disk. This causes the radial drift of the 
dust in the midplane to decrease due to collective effects. 
One would now intuitively think that a continuously de- 
creasing h/H ratio leads to smaller and smaller radial drift 
velocities but this is not the case for the following reason: 
When the ratio h/H decreases the vertical gradients of 
the azimuthal gas and dust velocities increase. Therefore, 
more angular momentum is transported in the higher re- 
gions of the disk. The midplane loses angular momentum 
which directly causes the radial velocity of the dust to 
increase. 

Finally both effects, the decrease in radial velocity due 
to collective effects and the increase of velocity due to the 
vertical transport of angular momentum, seem to cancel 
each other (or at most result in a negligible change in 
radial velocity) when small h/H ratios are considered. 



3.3. Effects of non-linear dynamics 



ratio e — Y.d/'Sg exceeds 0.40. For higher turbulence pa- 
rameters the critical radius interval decreases slower with 
higher dust-to-gas ratios. We find that for a — 10~^ the 
interval disappears for a dust-to-gas ratio of e = 0.75. We 
conclude that removing the gas may be a possibility to 
preserve mm to cm particles in the outer part of the disk. 

3.2. Turbulence parameter q 

Little attention was given to the turbulence parameter q 
until now (cf . Appendix lA.ip . A certain diffusion coefh- 
cient of the gas may be realized by big gas eddies moving 
slow or by small gas eddies moving fast. These two extreme 
cases are represented by g = 1 and q — 0, respectively. To 
illustrate how strongly q may influence the thickness of 
the dust layer we consider the following numerical exam- 
ple. We assume a Stokes number of unity and a turbulent 
a parameter of 10"'^. For the extreme case q — we cal- 
culate a dust scale height of h/H = 10~^ and for q — 1 
we obtain h/H = 3 x 10^^. These two dust scale heights 
differ by a factor of 30 which possibly influences the drift 
time scales. 

We calculated the effect of the turbulence parameter q 
on the drift velocity for two different turbulence parame- 
ters a, i.e. a = 10"^ and a = 10^^. The Stokes number 
was chosen in a way the the drift velocity for g = 1/2 
corresponds to a timescale of 1 Myrs. For a — 10^^ and 
a = 10^^ this implies St — 225 and St — 235, respectively. 
The results of these simulations are shown in Fig. [TH] 

This figure shows that the integrated drift velocities 
of the dust vary by 3% for a = 10~^ and by 1% for a — 
10^^ when q is changed from zero to unity. We find that 
for higher a parameters this variation is always less than 



So far we have considered different equilibrium states that 
would potentially allow solid particles to reside at large 
orbital radii for a longer time than a single test particle. 
It is however also possible for dynamical effects, such as 
spiral arms, turbulence or vortices, to reduce the radial 
drift. 

Dust particles are forced to climb up the local gas pres- 
sure gradient. In the simple case of a gas pressure that falls 
monotonically with radius the particles fall into the inner 
disc, but the particles may end up in a ny local gas oyerden 



sity that they encounter on their way (jKlahr fc Linl . 12001 



Haghighipour fc Bossllioosh . If the disc is massive enough 
to be gravitationally unstab le, its spira. 1 arms may act as 



such local density maxima (jRice et al.l . 120041 ) . Transient 



overdensities that occur in magnetorotational turbulence, 
i n a way very elongate d vortices, have the same effect 



(jJohansen et al.l l2006b[ ). slowing down radial drift by a 



factor of two. The important parameter for reducing the 
overall radial drift is the life-time of the gas overdensity. 
Spiral arms would from this perspective be a good candi- 
date, since turbulent overdensities tend to live no longer 
than a few local orbits of the disc. 

The coupled flow of gas and dust is in itself un 



stable to the streaming instability (jYoudin fc Goodmanl . 
2005) , lea ding to particl e clumping in the non-linear state 
l2006bh . These local dust overdensities 



(Johansen et al 



drag the gas along with their orbital motion, thus reduc- 
ing the sub-Keplerian head wind and the radial drift. The 
effect of the streaming instability on t he radial drift can be 
as hi gh as a factor two in reduction (jJohansen fc Youdinl . 
2007h . 
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3.4. Temperature and surface density profiles 

We also investigated to which extent the radial drift time 
scales of the dust particles depend on the temperature and 
the surface density profile. These quantities were found to 
play a minor role which can be reasoned as follows. If 
the temperature is decreased by a factor of 2 then the 
maximum radial drift velocity wn decreases by the same 
factor (cf. Eq. IA.22P . This means that if the temperature 
at 100 AU is decreased from 20 K to a rather low value 
of 10 K then the drift time scales would only increase by 
a factor of 2. Since the radial drift problem spans at least 
one order of magnitude a change in temperature does not 
provide a solution. A change in the surface density index 
from (5 = 0.8 towards d ~ changes the maximum drift 
velocity wn by a factor of 1.5 according to Eq. (|A.22[) which 
also does not solve the problem. 

4. Summary 

We discussed two different subjects. First, we investigated 
the velocity structure of the dust layer and the gas layer. 
We calculated typical radial drift velocities of the dust in- 
cluding several effects, i.e. collective effects and effects due 
to turbulent viscosity. The results of this first issue enabled 
us to estimate characteristic drift time scales of the dust 
which represents the second area of interest. These calcu- 
lations were performed in order to explain interferometric 
millimeter observations that indicate the presence of mm 
to cm size particles in the outer part of disks with ages up 
to 10 Myrs. In the following we will briefly summarize the 
results. 



conclude that collective effects play the major role in the 
parameter regimes considered in this paper. 

It turned out to be neccessary to include radial vis- 
cosity terms of gas and dust. If these terms are neglected, 
the numerical solutions showed very steep and unrealistic 
velocity gradients. Radial viscosity, once included in the 
considerations, reduced these velocity differences which 
considerably changed the radial velocities of gas and dust. 

For turbulence parameters a larger than 10"** collec- 
tive effects as well as effects of turbulent viscosity play a 
minor role. The dust is stirred up away from the midplane 
in such a way that the dust-to-gas ratio is too low for these 
effects to play a role. The radial drift velocities of the dust 
in this case were reduced by a factor of 1.3 at most. For 
this reason, both effects can be disregarded in this regime. 

We compared the results of our numerical solutions 
to the predictions of the so-called 'plate drag approxi- 
mation'. It was shown that both models roughly agree 
if the Stokes number of the particle is about unity. For 
Stokes numbers higher (lower) than unity the two mod- 
els predict different radial velocities. Previous work al- 
ready put the 'plate dra g approximation' into question 
(jYoudin fc Chiang . 120041 ) with similar conclusions. 

Apart from radial drift behaviour some structural as- 
pects of the gas layer and dust layer were also investigated. 
We calculated the extent of the vertical region in which gas 
and dust affect each other. This calculation showed that 
the gas layer behaves like a classical Ekman layer if Stokes 
numbers larger than unity are considered. For well coupled 
particles, i.e. for particles with Stokes number smaller than 
unity, the two descriptions differ significantly. 



4.1. Velocity structure of gas and dust layer 

First we discussed the individual parti cle drift of the dust 
(jWhippld . [l973 IWeidenschillintd . Il977l ) . In this model the 
dust particles are considered to be independent and the 
gas is assumed not to be affected by the dust at all. This 
scenario involves no vertical structure of the dust since 
only individual dust particles are considered. 

We discussed the influence of collective behaviour on 
the radial drift velocities of the dust. Collective effects 
take place when the dust-to-gas ratio exceeds unity, e.g. 
when the dust sediments into a thin midplane layer due 
to low turbulence. By assuming an initial dust-to-gas ra- 
tio of 0.01 we found that the local radial drift velocity of 
the dust in the midplane of the disk may be reduced by 
a factor of 100 for small turbulence parameters and cer- 
tain particle characterictics. The vertically averaged radial 
dust velocity is at least an order of magnitude lower than 
the corresponding individual particle drift. 

The additional inclusion of turbulent viscosity in our 
considerations increased the radial drift of the dust due to 
the vertical transport of angular momentum. This increase 
in radial velocity was not as significant as the decrease by 
including collective effects. Turbulent viscosity enhanced 
the radial drift by a factor of 2 at most. Therefore, we 



4.2. Radial drift time scales 

We investigated the radial drift of dust particles with re- 
spect to collective behaviour and effects of viscosity un- 
der standart conditions, i.e. an initial dust-to-gas ratio 
of 10~^, standard turbulence parameter q, etc. We found 
that a possible solution for the drift problem are very 
high/low disk masses. Moreover, high porosity particles 
(/ < 0.1) in low turbulent disks with masses > 0.2 AI^, 
provided drift time scales > 2 Myrs. 

However, except for these two cases the drift time 
scales of the dust turned out to be far too short to ex- 
plain the millimeter observations. For this reason we in- 
vestigated several possibilities in order to increase the drift 
time scales. We found that removing the gas from the disk 
may increase the time scales up to more than 2 Myrs. For 
example, for a turbulence parameter of a = 10^^ we found 
that cm size particles remain in the outer parts of the disk 
for more than 2 Myrs if disk masses < 0.05Af^ or > 0.2^^, 
are considered and only 5% of the gas is left. If the disk 
is more turbulent more gas has to be removed in order to 
keep the particles in the outer parts for longer periods of 
time. 

Gas might be removed in the later stages of disk evo- 
lution when photoevaporation sets in. While the gas is 



16 



Brauer, DuUemond, Johansen, Henning, Klahr and Natta: Radial drift of mm-cm size grains 



evaporated from the disk the dus t part icles > 20 fim 



(|Takeuchi fc Linl . l2005t iBallv et all |2005[ ) remain in the 



outer parts of the disk. However, current photoevapora- 
tion models remove the gas from the disk rather abruptly 
( Alexander et al. . 2006[ ). This would lead to rather high 
relative particle velocities in the disk. Hence, the dust par- 
ticles would fragment and destructive collisions would play 
an important role. The centimeter particles would be de- 
stroyed in short time scales which is not in agreement with 
the observations. 

We also investigated the effect of the turbulence pa- 
rameter q on the drift time scales and we found that this 
parameter plays a minor role. Varying this parameter be- 
tween zero and unity, the most extreme cases possible, 
changed the drift velocities by 3% at most. 

We did not consider the growth of the particles in 
this paper. The growth time scale of the dust to reach 
a centimeter in size could be of the order of 1 Myr con- 
sidering the outer parts of the disk. A model including 
both processes, radial drift and coagulation, would clar- 
ify this issue. Recent work about the drift time scales 
in comparison to the growt h time scales was done by 
iKlahr fc Bodenheimerl According to their calcu- 

lations, the dust would long have drifted away from 100 
AU before the particles even reach the size of about a 
centimeter. 

Moreover, the effect of fragmentation could change the 
situation. Even though a process that destroys particles is 
included it could finally speed up the process of coagula- 
tion: Fragmentation leads to a permanent amount of small 
particles as a result of collisions. These small particles may 
be swept up by larger particles due to their high relative 
velocity. Although some particles are destroyed, the final 
sum of both effects might lead to an accelerated growth. 
These effects, the radial drift, the dust particle coagula- 
tion and the effect of fragmentation, will be the topic of a 
forthcoming paper. 

4.3. Observational implications 

We have shown that thin midplane layers in low turbulent 
disks can conceivably explain the presence of the observed 
large amounts of mm/cm sized grains if a significant frac- 
tion of the gas is removed, for instance through photoe- 
vaporation of the gas. Therefore, we should also discuss 
whether the presence of such thin midplane layers is con- 
sistent with observations of e.g. edge-on disks. 

The infrared emission from protoplanetary disks orig- 
inates from smaller 3 mm) dust grains. These grains 
must be smaller than 3 micrometer, as can be inferred 
from the presence of a 10 micron silicate feature in emis- 
sion in most sources. Even with relatively little verti- 
cal mixing (low turbulence) the very smallest dust par- 
ticles can still be mixed up to intermediate he ight above 



II in the Meeus et al.l (|2001r ) classification. Interestingly, 
Acke et al. ( 2004) have shown that there is indeed a cor- 
relation between the presence of large grains in the outer 
regions of disks around Herbig stars and the type of SED 
of the disk: disks with large amounts of large (mm/cm) 
grains appear on average to also have SEDs consistent 
with flatter disk geometry. This seems to substantiate at 
least qualitatively the idea of low turbulent disks. 

For larger grains, which can be observed at mm/cm 
wavelengths, there have not been observations of edge-on 
disk with very thin mm disks so far. However, this has two 
reasons. The first reason is that if the mm/cm disk is very 
thin, then the chance to observe it sufficiently precisely 
edge-on is very small. This reduces the number of poten- 
tial candidates for such measurements drastically. The sec- 
ond reason is that current state-of-the-art interferometers 
do not yet have the spatial resolution to resolve such thin 
disks. For instance, the Butterfly Nebula (a well-studied 
ne arly perfectly edge -on disk), was resolved with OVRO 
bv IWolf et all (|2003h . but the vertical extent of the ob- 
served disk was as large as the beam size, i.e. unresolved 
in vertical direction. 



5. Conclusion 

We find that the radial drift time scales for millimeter- 
sized grains observed in the outer (~ 100 AU) regions of 
protoplanetary disks are of the order of 10^ years, whereas 
the age of these disks is generally between 10^ and 10^ 
years. So, according to theory, the particles that are ob- 
served should have long ago drifted into the star, which 
is evidently inconsistent with observations. We have at- 
tempted to resolve this mystery by investigating what 
happens in the case of very low-turbulence disks in which 
the dust gathers in a thin midplane layer. The hope was 
that collective effects of the dust (by coupling to the gas) 
could slow down the drift. We find that this reduction 
factor is sufficient if either very high/low disk masses are 
considered or the particles have high porosities. In any 
other case this reduction factor is by far insufficient to 
resolve the problem. Moreover, we find that by removing 
a substantial fraction of the gas from the disk the drift 
time scales can be sufficiently enhanced. Another possi- 
bility is particle t r appiii g in local gas pressure maxima 
([Klahr fc Henning . 119971 ) , but this topic was beyond the 
semi-analytical model of this paper. 



the midplane (jDullemond fc Dominikl . l2004al ). although 
we admit that the disk should look significantly less 
flared in such a case, i.e. the disk should be of Group 
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Appendix A: Theoretical background 

The main ingredients of the models arc the stellar param- 
eters, the structure and mass of the disk, the description 
of the turbulence, and the description of the interaction 
between the dust and the gas. In principle, for the ra- 
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dial drift velocity we do not need knowledge of the tur- 
bulence (except for th e turbulence effects described by 
Johansen et al. I (|2006bh which might reduce the drift by 
a factor of up to 2). However, turbulence sets the thick- 
ness of the dust midplane layer, and hence the density of 
the dust, which affects the drift velocity once this den- 
sity comes close to, or exceeds the gas density. If there 
would be no turbulence at all, the dust particles would 
settle to the midplane and form a razor-thin layer qual- 
itatively like the rings of Saturn. The more turbulence 
there is the harder it is for the dust to form a thin layer 
since it is mixed up and again transported in the higher 
regions of the disk. Turbulence is not yet understood in 
detail and there is var ious literat ure about this challeng- 
ing topic (IVolk et al.l. IT 980;_ Schra pler fc Hennhlgl . 12004 : 



Balbus fc HawlevLll99lL " 



9981). In this paper we will make 



use of the a-prescription ( Shakura fc Sunvaev , 1973f ) . This 
specific way of regarding turbulence is somewhat superfi- 
cial since it ignores most details of turbulence. However, 
the advantage is that it makes turbulence manageable 
without extensive hydrodynamical simulations. 

A.l. a-prescription for turbulent gas 

In this section we will give a short overview about our 
a-prescription of particles in a turbulent disk. Turbulence 
mixes things up and therefore acts like a kind of diffusion. 
The diffusion coefficienl|f| D can be written as a product 
of a velocity scale Vq and a length scale Lq, 



D^VnXLn 



(A.l) 



A typical velocity scale is the isothermal soundspeed Cg = 
fcTgas//^, with Tgas the gas temperature, k the Boltzmann 
constant, fi the mean molecular weight of the gas (which 
we take 2.3 times the proton mass for a standard mixture 
of molecular hydrogen and atomic helium) . A characteris- 
tic length scale is given by the pressure scale height of the 
gas H = Cs/VLk, where fifc = ^jGM^/r'^. Regarding these 
typical scales at a certain radius we can alternatively ex- 
press the diffusion coefficient for the gas as follows 



D = ac,H 



(A.2) 



The value of the parameter a reflects the amount 
of turbulence in the disk and it ranges from 10~^ 
( Weidenschillingl . Il980( ) up 



to 



r anges 
10^^ ( Hartmann et al 



19981 ). Extensive hydrodynar nical simulations sh o w typ - 



ical a values of about 1 ^ ()Brandenburg et al.l (|l995f ). 
Johansen fc Klahil (120051) ). In general we have to distin- 
guish between the a introduced to determine the vertical 
distribution of the dust and the Uacc linked with the ac- 
cretion process. Both quantities are related by aacc =PtQ; 
where Pt denotes the Schmidt n umber . We will assume 

1975 



Pt= 1 in this paper ( 


Launder 


Johansen fc Klahr, 2005 


). 



McCombl . Il990l: 



* All symbols used in this paper are listed in a table at the 
very end. 



There is an ambiguity with regard to this formulation 
till now since a does not provide any information about Vq 
and io- This can be seen in a better way by introducing 
a turbulence parameter q 



D = aCsH 



a^Cs X a 



(A.3) 



which ranges between and 1. Now we can identify 
the terms of equation (|A.1|) by defining Vq — a'^Cg and 
Lq = a^~'^H. The physical interpretation of q is as fol- 
lows. A certain diffusion can be realized by big eddies 
which move very slow (q = 1) or by small eddies mov- 
ing very fast {q = 0). Both possibihties result in the 
same diffusion-coefficient. Norma lly q is set to be 1/2 
(e.g. Schrapler fc Henningl ( 20041 )) but there are other 
values possible. For example, in self-induced turbulenc e 
q tends to be smaller than 1/2 ( Weidenschillind . 2006f ). 
Another possibility are big convective cells of scale H 
(|Klahr fc Hennind . Il997[) which would imply q = 1. The 
fundamental importance of knowledge about q can be un- 
derstood by calculating the turn-over-eddy-time Tod 



Ted 



Lo 
Vo 



27r 



vl-29. 



(A.4) 



where Ucd is the eddy turn-over frequency. Comparing the 
two extreme cases q — and q — I the timescale Tod 
changes by a factor of . Taking a typical value of a = 
10"'^ means that these time scales differ by 6 orders of 
magnitude. In the next section we will see that the effect 
of the turbulent gas on the dust, and therefore on the dust 
scale height, is highly dependent on this timescale Ted- For 
this reason knowledge about q is an essential requirement. 

A.2. Viscosity of gas and dust 

In this paper we will assume that the value of the tur- 
bulent viscosity of the gas is basically the turbulent dif- 
fusion coe fficient of the gas ea. (IA.2I) . which corresponds 
to Pt = 1. Johansen fc Klahil (120051 ) found that the ratio 
between these two quantities ranges from 0.8 up to 1.6 in 
MRI turbulence, depending on the direction, i.e. vertical 
or radial. 

We will also assume a certain turbulent viscosity of 
the dust. This quantity can be expressed as the product 
of a characteristic mixing length and a relative turbulent 
velocity of the particles 



(A.5) 



For Stokes numbers greater than unity the relative tur- 
bulent veloci ty between the dust particles is g i ven b y 



= Vo/St (|Volk et all Il980l : ISuttner fc YOTkl l200ll) . 



The velocity Vq denotes the turbulent velocity of the 
largest gas eddy introduced in the ffi'st sections. With the 
mixing length Lmix — TgVt the viscosity of the dust is given 

by 

" for St > 1 , (A.6) 



D 

St 



where we used the Eq. (iX2|) and [KJO\i . 



Brauer, DuUemond, Johansen, Henning, Klahr and Natta: Radial drift of mm-cm size grains 



19 



Now, for Stokes numbers smaller than unity the dust 
particles are well coupled to the gas, i.e. both disk com- 
ponents, gas and dust, behave more like one single fluid 
than two different types of matter. We have already seen 
that the diffusion coefficient for the dust in eq. (|A.14|1 
matches the diffusion coefficient for the gas in the St < 1 
regime. We will assume the same behaviour with regard 
to the viscosity, e.g. that the dust viscosity equals the gas 
viscosity for small St. Considering this, the dust viscos- 
ity in both regimes, S t > 1 and St < 1, is then given by 
( Schrapler fc Henningl . [2004 



formula: 



D 



Vd = 



1 + St 



(A.7) 



A. 3. Self-induced turbulence 



In this paper we will vary a between a rather high value 
of 0.01 down to a rather low value of 10~^. As has been 
discussed by many authors (and will also be discussed be- 
low) we expect the following behaviour: For low a the 
dust sediments into a very thin midplane layer. However, 
there is a limit on how thin this layer can become, or 
in oth er words: how low the level of turbulence can be- 
come. IWeidenschilliMi (|l979l ) has shown that a shear- 
instability between the dust layer and the gas induces a 
weak, but non-negligible level of turbulence. This is called 
'self-induced turbulence'. 

In principle this self-induced turbulence can be de- 
scribed by the (a,g)-formalism discussed above, as long 
as we allow that both a and q depend on height above the 
midplane z and radial distance to the star r. Determining 
the a and q of this self-induced turbulence is a complex 
matter, and although enormous progress has been made 
in this field, until now there are still many open issues 
( Cuzzi et al. . 1993t Weidenschilline . 20061 : Johansen et al 
2n06ah . 



In this paper we will not explicitly address the issues 
of self-induced turbulence: we will keep a and q as param- 
eters of the model (assuming that they do not depend on 
z or St). However, in appendix[B]we will roughly estimate 
the level of self-induced turbulence to obtain a lower limit 
to the a that we are allowed to use. 

A. 4. Stokes number 

A moving particle in a gas at rest loses a significant frac- 
tion of its momentum within a time called stopping time 
Ts ■ This time depends on the friction between the particle 
and the gas. A strong friction means a small r^, and vice 
versa. The friction depends on the particle cross section 
(Tp = TTo^ and, therefore, the particle radiu^ a, the rela- 
tive velocity Vp with respect to the gas and the properties 
of the gas (mainly gas density pg, isothermal sound speed 
Cs and molecular mean free path I). 

For particles, with a size smaller than the molecular 
mean free path, the friction can be expressed by a simple 



(A.; 



This is the "Epstein drag law" . In this regime the stopping 
time equals: 

_ nipVp Ep. psa 



PgCs 



(A.9) 



where nip is the particle mass, which can be expressed 
with the particle material density ps as = (47r/3)psa'^- 

For particles larger than the mean free path the drag 
law is much more complex. This regime is characterized by 
the "Stokes drag law" . In this paper we focus on particles 
that are always smaller than the mean free path and can 
ignore the Stokes regime. 

If the stopping time Ts is much smaller than the turn- 
over-eddy time Ted, the particles are strongly coupled to 
the gas having the same motions and the same behaviour 
with regard to diffusion. When exceeds Tod by far, the 
dust decouples from the gas and is hardly influenced by 
the turbulence of the gas. 

The stopping time characterizes the dynamic proper- 
ties of the particle as it moves through the disk. Therefore, 
we can replace all microphysical particle properties like 
a, Ps, nip by Ts- Particles with vastly different properties 
(e.g. size), but the same Ts behave entirely the same. 

Instead of using the stopping time Tg, an even more 
convenient parameter is the so-called "Stokes number" 
Sti. It is defined by: 



Str 



Ted 



Ts^ka'"-^ 



(A.IO) 



The particles are strongly coupled for St^ <C 1 and hardly 
influenced by the gas for St/, ^ 1. For the case q = 1/2 
the Stokes number does not depend on a. We will now 
introduce the Stokes number St by St=StL(9 = 1/2) = 
Ts^k and we will formulate all the particle properties in 
terms of St and St l ■ 

A. 5. Thickness of the dust layer 

Let us now turn to the calculation of the thickness of the 
dust layer. This thickness is determined by an equilib- 
rium between dust which settles towards t he midplane and 
diffusion which stirrs the dust up again (IDubrulle et al. , 
ll995 ljSchrapler k Hennind . 120041: iDullemond fc Dominikl . 
i2004br ). The settling can be described by the equation of 
motion for a dust particle, 



-fi^z z . 

Ts 



(A.ll) 



For St > 1/2 the particle describes a damped oscillation 
around the midplane, corresponding to an orbit inclined 
with respect to the midplane. For St ^ 1/2 the particle 
is so well bound to the gas, that the downward motion 
equals to an equilibrium settling motion with magnitude: 



We will always assume the particles to be spherical. 



Wsott 



— — zStrik 



(A.12) 



20 



Brauer, DuUemond, Johansen, Henning, Klahr and Natta: Radial drift of mm-cm size grains 



SO that the second order differential equation Eq. (|A.lip 
reduces to a first order differential equation: 

i = -Wsott • (A.f3) 
The diffusion of the dust is characterized by 

Dh = — ^ , (A. 14) 



l + Sti 



which is the diffusion coefficient of the gas corrected by a 
factor including the coupling between dust and gas. This 
expression decreases with higher Stokes number since the 
dust is not longer affected by the gas. From the numbers 
Dd and Tgott, representing settling and diffusion, we can 
construct a length scale by ft.^ = DdTsctt- This leads to the 
expression for the dust layer thickness: 



= Dd max(Tsctt, l/2f^k) 



(A.15) 



We resticted the settling time scale to be at least half 
an orbital time scale. For St larger than unity the mo- 
tion of a dust particle above the midplane corresponds to 
a damped inclined orbit. The velocity towards the mid- 
plane of the disk can not be higher than the projected 
Kepler velocity V\^z/r and, hence, the time for the par- 
ticle to cross the midplane not significantly smaller than 
an orbital time scale. However, the dynamics of bodies in 
quasi-Keplerian orbits in turbulence may be not well de- 
scribed by diffusion, and the authors are aware that the 
approach Eq. (|A.15p should be used with caution. Also 
the vertical distribution of the dust for the case of in- 
clined orbits (i.e. for large Stokes number) is not Gaussian 
but rather bimodal, but we will ignore this effect. With 
Eqs. (jA.2[IA.l"4|) . as well as the expression for the gas scale 
height H = Cg/rik, one can rewrite Eq. (|A.15|) as: 



min(St,l/2) (1-1- Sti) 



(A.16) 



This is the most general description of the dust layer. 
The Stokes numbers St and St^, allow us to calculate the 
thickness of the dust layer and, therefore, the dust mass 
densities. With these densities we are able to calculate 
radial drift velocities to estimate particle radial drift times 
in the outer parts of the disk. 

A. 6. Gas and dust density 

With the results of the last section we now present the 
mass densities of gas and dust. We assume the disk to 
be isothermal in the z-direction and we use a thin disk 
approximation: z <C The vertical distribution of the gas 
is then given by 

2i/2 



2nH 



exp 



(A.17) 



The scale height of the gas H is given by H — Cg/r^k- The 
density distribution of the dust is 

Pd = ^exp('-l^V (A.18) 



2TTh 



wherein h denotes the scale height of the dust we esti- 
mated in the last section and Ed ist given by eoSg. The 
quantity eo in this expression is the initial dust-to-gas 
mass ratio and it is usually set to 0.01. Further investiga- 
tions will show that the drift time scales depend strongly 
on this parameter which means that eq is one of the main 
quantities describing our system. Let us assume that the 
surface density Eg of the gas follows a power law 



Eg(r) = Eo 



and the temperature T does so as well: 



T(r) 



(A.19) 



(A.20) 



We choose the power law index 6 to be 0.8 following 
iKitamura et al. (|2002 ). Moreover, we assume the temper- 
ature To to be 200 K, corresponding to a passive disk ir- 
radiated under an angle of 0.05 degrees around a T Tauri 
star with a surface temperature of 4000 K and i?* — 2.5 
Rq. The power law index of the temperature ^ is set to 
be 1/2. The mass of the disk is given by 



disk 



Yig{r)rdrdip 



27rEor| 
2-6 



2-5 _ 2-Sl 
max cv I 



(A.21) 

which we take as a parameter of our models. The size 
of the disk Tmax is set to be 150 AU following the 
observati onal r e sults of the Taurus- Auriga region by 
Rodman n et al. ( 2006h . The dust evaporation radius rev 
is 0.03 AU. Throughout the paper we will give the mass 
of the disk in terms of the central stellar mass Mi, which 
we assume to be 1 Mq. 

If we consider Herbig Ae/Be stars instead of T Tauri 
stars then the disk mass and the outer disk radius in- 
creases (Natta, 2004). Herbig stars are usually a factor of 
4 heavier and have outer disk radii of 300 AU and more. 
However, Eq. (jA.21|) shows that the surface density Eq 
is hardly affected by this change. The disk is more mas- 
sive, but since the disk extent increases as well the actual 
amount of mass per unit area in the disk does not change. 

In this paper we use a standard velocity wn- This quan- 
tity is the velocity by which the gas moves azimuthally 



slower than Kcplcrian veloc ity Vk , i.e 



is given by (Weidenschilling. Il977 : Nakagawa et al.l . ll98(ih 



V\<-VN. It 



UN 



drPg 

'2pgrik 



2Uk 



2 + 2+^ 



(A.22) 



Appendix B: Self-induced turbulence 

There is a limit of how thin the dust layer can be, i.e. how 
low the level of turbulence can become. The calculations in 
this paper are physically realistic as long as the turbulence 
parameter a does not drop below this value. In this part 
of the appendix we present the lowest possible amount of 
turbulence in terms of the turbulence parameterisation a 
and q. 
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The reason for a minimal amount of turbulence is the 
following: when the dust settles towards the midplanc and 
the dust-to-gas ratio exceeds unity the gas is carried along 
with the dust. Since the gas above the midplane still moves 
somewhat sub-keplerian the resulting shear can destabi- 
lize the layer causing turbulence. This instability is called 
Kelvin-Helmholtz instability. The turbulence in turn again 
diffuses the particles until a steady state is reached. 

In the next two subsections we will translate this 
effect into the a-prescription. To do this, we have to 
consider two different cases. The Stokes number Stc 
which divi des these two cases is roughly given by 10^^ 
dCu^zi fc Weidenschininel . l2006h . Let us first focus on the 
St < Str case. 



be calculated to be 



1 



q = 



1 - loga 



2tt 
Ro 



(B.4) 



Taking Ro=25 we get q — 0.45 which is quite close 
to the value 0.5. If we assume a /? of about 10~^ at 
100 AU then the self-induced turbulent a parameter 
in this case is roughly given by 10~^. According to 
ICuzzi fc Weidenschilline ( 2006h . the value of the Rossby 
number is uncertain up to a factor of 2-3. This means that 
the a-value could be even lower and also in the order of 
~ 10~^ at 100 AU in the disk. This might be also the case 
if higher stellar masses or lower temperatures are consid- 
ered. 



B.l. St< Stc 

In this regime of small particles the thickness of the layer 
is almost inde pendent of particle properties and given by 
(|Sekival . [l998h 



h r^H 



(B.l) 



The quantity Ri is the critical Richardson number. This 
number, which is usually assume d to be constant , is ap - 



proximately given by Ri = 0.25 (jChandrasekharl . 119611 ). 
Recent publications call this assumption in question 
( Gomez fc Ostrikeiilioosh and show that the Richardson- 
number required for marginal instability is eventually 
higher than the traditional value of 0.25 depending on the 
Stokes number. However, since the results of our simula- 
tions turned out to be hardly dependend on Ri this effect 
may be neglected in the context of radial drift velocities. A 
comparison between (|B.1|) and (|A.16|) yields the turbulent 
a parameter for this case. 



a = /?RiSt(l + StL) 



(B.2) 



The quantity [3 is simply the squared ratio between the 
isothermal sound speed Cg and the Kepler velocity 14- At 
100 AU the quantity (3 is roughly of the order of 10~^. 
This leads to a turbulent a value which is at most of the 
order of a ~ 10^^ assuming q — 1/2 and Ri — 1/4. 



B.2. St > Stc 

In this second regime the diffusion coefficient in self- 
induced turbul ence for the gas is g i ven by D = 
(^yk)Vf^kRo^ (|Cuzzi fc Weidenschillingl . l200(ih . A direct 
comparison between the last equation and Eq. (|A.2p shows 
that the turbulence parameter is a constant: 



Ro^ 



(B.3) 



To estimate the q parameter we have to focus on the so- 
called Rossby number Ro that represents the ratio be- 
tween the turn-over-frequency of the largest eddie Wed and 
the Kepler-frequency, so that u^d = rifcRo. With this re- 
lation and equation (IA.4|) the turbulence parameter q can 
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Variable Explanation 

a radius of the dust particle 

rup mass of the dust particle 

(7p geometrical cross section of the dust particle 

Vp velocity of the dust particle 

/ filling factor of the dust particle 

r distance to the central star from a point in the midplane 

z height above the midplane 

Cs isothermal soundspood 

fik, Vk Kepler frequency, Kepler velocity 

H = Cs/f2k gas scale height 

h dust scale height 

Eg , Ed surface density of the gas and the dust 

Pg, pd gas and dust density 

2 

Pg = CgPg gas pressure 

eo, e initial and current dust-to-gas ratio 

V' = 1/(1 + e) reduction parameter in the NSH solution 

77 = — 9rPg/2pgrf]k ratio between radial pressure force and gravitational force 

Tev evaporation radius 

Tmax outer radius of the disk 

mass of the central star 

Mdisk mass of the disk 

vn = ??Vk typical radial drift scale 

a, q turbulence parameter 

Ts stopping time of the dust particle 

St = TsfJk Stokes number of the particle 

St_L = StQ^''"^ modified Stokes number 
Lo, Vb, Wed, Ted length, velocity, frequency scale and time of the turbulent gas motions 

£) = i/g, Di = I'd diffusion coefHcients/viscosity of gas and dust 

Re* critical Reynolds number 

Wr, Wip radial and azimuthal velocity in the single particle drift equations 

u"^^^ radial velocity of the dust in the NSH solution 

u^^^ vertically integrated radial velocity of the dust in the NSH solution 

radial and azimuthal velocity of the gas (solution of the NSe) 

Wr, Wip radial and azimuthal velocity of the dust (solution of the NSe) 
vertically integrated radial velocity of the dust (solution of the NSe) 

A vertical lenght scale in which dust and gas effect each other (NSe) 

Ag k'ligiit Kcale of the elassieal Eknian layer 

(5p drift iiithiccd l)y tlio iuclusiou of viscosity 



Table B.l. Variables used in the course of this paper. 



